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Abstract —Sustained oscillations observed in power systems can 
damage equipment, degrade the power quality and increase the 
risks of cascading blackouts. There are several mechanisms that 
can give rise to oscillations, each requiring different countermea¬ 
sure to suppress or eliminate the oscillation. This work develops 
mathematical framework for analysis of sustained oscillations 
and identifies statistical signatures of each mechanism, based on 
which a novel oscillation diagnosis method is developed via real¬ 
time processing of phasor measurement units (PMUs) data. Case 
studies show that the proposed method can accurately Identify 
the exact mechanism for sustained oscillation, and meanwhile 
provide insightful information to locate the oscillation sources. 

Index Terms —Power system stability, oscillation diagnostics, 
phasor measurement units, weakly damped oscillation, limit 
cycle, Hopf bifurcation, forced oscillation. 

I. Introduction 

USTAINED low frequency oscillations are one of major 
concerns to power system operation. Oscillations cause 
problems for power quality and can potentially damage power 
grid equipment or activate protective equipment. In most 
severe scenarios, growing oscillations may lead to catastrophic 
blackouts [1]. 

Generally speaking, there are three mechanisms of power 
system oscillations. Firstly, the oscillation can appear due 
to weak damping arising from high-gain fast exciters, long 
transmission lines, or high transmission powers [2]- [4]. This 
kind of oscillation can be quenched by appropriate power 
system stabilizer (PSS) tuning, intertie line controls, and 
generator power reduction. Second mechanism of oscillation 
appearance is attributed to supercritical Hopf bifurcation such 
that a stable limit cycle is born [5]- [10]. Unlike the weakly 
damped oscillation which can be analyzed by linearizing the 
system, the limit cycle is an essentially nonlinear behavior 
which only exists in nonlinear systems. The emergence of 
limit cycle implies that the stability property of the system 
has changed, and in order to push the system back to stability, 
operation point may need to be reset. In addition, another 
mechanism of oscillation is forced oscillation, which is excited 
by external periodic disturbance including cyclic loads, control 
loops in power plants, diesel generators, etc. [11]- [15]. 
The most effective countermeasure is to locate and separate 
the external disturbance from the system. In summary, there 
are different mechanisms of power system oscillations, and 
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different countermeasures need to be adopted accordingly. 
Hence, detecting power system oscillations and diagnosing the 
corresponding mechanism are of utmost importance in power 
system monitoring, and of great significance to ensure the 
secure operation of any power system. 

Phasor measurement units (PMUs) have been widely de¬ 
ployed in power grids to provide system states and dynamics 
in real time [16] [17]. The high quality PMU data provides in¬ 
valuable information to enable oscillation diagnosis. Regarding 
the weakly damped oscillation, different methods to identify 
oscillation modes and damping ratios are proposed including 
Prony analysis [18], frequency domain decomposition analysis 
[19], subspace identification method [20], robust RLS methods 
[21] and so forth. Hopf bifurcation has been studied in terms 
of dynamic stability. Hopf bifurcation not only involves with 
the oscillation, but also closely relates to the voltage stability. 
The previous studies show that voltage collapse may arise 
from the existence of Hopf bifurcation, which is prior to the 
appearance of saddle-node bifurcation [5]- [10]. In addition, 
forced oscillation has also been studied by exploring frequency 
domain techniques [11]- [15]. It has been shown that if the 
forced oscillation is close to the natural frequencies, resonance 
may be observed which leads to severe consequences [3] 
[14] [15]. Even though each type of the oscillation has been 
widely investigated, little effort has been made to distinguish 
the three mechanisms from time series data such that the 
right control actions can be executed. The challenges of 
oscillation diagnosis include low signal-to-noise ratio (SNR) 
especially when the amplitude of oscillation is small while 
load fluctuation and noise intensity are large, as well as similar 
characteristics of time series data. Specifically, the time series 
data of oscillations with different mechanisms all looks alike 
to each other. 

This paper proposes a novel method to diagnose the mech¬ 
anisms of sustained oscillations using PMU data. A unified 
mathematical framework to describe sustained oscillation is 
developed, under which the statistical signatures of different 
models are explored. It is shown that sustained oscillations 
of different mechanisms exhibit distinct statistical signatures 
including the kurtosis and the power spectral density, based on 
which an oscillation diagnosis method is developed. Numerical 
examples show that the proposed diagnosis method is able to 
accurately identify the exact mechanism for oscillation even 
when the SNR is low, and meanwhile provides insightful 
information to locate the oscillation sources, which are crucial 
to implement the right control actions in a timely manner. 
Note that the goals and focus of the proposed method and 
conventional waveform analysis techniques are different and 
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complementary to each other. The proposed method is to 
diagnose different oscillation mechanisms and locate the po¬ 
tential sources, whereas the conventional waveform analysis 
techniques focus on analyzing damping, frequency, oscillation 
mode shapes, etc. By nature, the traditional waveform analysis 
does not distinguish between different oscillation mechanisms. 

The paper is organized as follows. Section II describes a 
unified framework under which the mathematical model for 
each oscillation mechanism is developed. Section III presents 
the statistical signatures for each oscillation mechanism based 
on which an oscillation diagnosis method is proposed. Case 
studies are presented in Section IV to demonstrate the accuracy 
and feasibility of the proposed method. Conclusions and 
perspectives are given in Section VI. 


II. MECHANISMS OF SUSTAINED OSCILLATIONS 
The power system dynamic model can be described as: 

i = f{x,y) (1) 

0 = g{x,y,u) (2) 

Equation (1) describes dynamics of generators and their asso¬ 
ciated control as well as load dynamics, and (2) describes the 
electrical transmission system and the internal static behav¬ 
iors of passive devices, f and g are continuous functions; 
vectors x € R"” and y £ R"« are the corresponding 
state variables (generator rotor angles, rotor speeds, etc) and 
algebraic variables (bus voltages, bus angles, etc), and u is 
the vector describing stochastic behaviors in real-world power 
systems. The stochastic perturbations can be originated by 
loads variations, renewable energy power injections, transient 
rotor vibrations of synchronous machines, measurement errors 
of control devices, etc. [22] [23]. In this paper, we are inter¬ 
ested in the stochastic perturbations like load variations and 
renewable energy power injections, which can be modeled as 
the Ornstein-Uhlenbeck process, which is stationary, Gaussian 
and Markovian [23] [24]: 

ii = —Cu + ai (3) 

where C = diagjai, a 2 ,..., Qf„„} which relates to the corre¬ 
lation times of the stochastic processes [23], is a vector of 
independent standard Gaussian random variables , is the 
intensity of noise. 

Linearizing (1) and (3), and also eliminating algebraic 
variables y from (2), the stochastic power system model takes 
the form: 


i) = Av -f aB^ 


(4) 


where v = \5x,5uY^, A 


fx fydy 9x fy9y 9u 

0 -c 


In the rest of the paper, we focus on the stochastic model 
described in (4), where u is a vector Ornstein-Uhlenbeck 
process. 

As discussed in the introduction, sustained oscillations of 
power systems can be produced by several vastly different 
mechanisms. In the sections below we discuss possible origins 
of oscillations and develop the corresponding mathematical 
models. 


A. Weakly Damped Oscillation 

Weakly damped oscillation can be classified into local 
mode oscillations, and interarea mode oscillations. Local mode 
oscillations are usually caused by automatic voltage regulators 
(AVRs) operating at high output and feeding into weak trans¬ 
mission networks; Interarea mode are associated with weak 
transmission links and heavy power transfers. Power system 
stabilizers (PSSs) are the most common means of enhancing 
the damping and suppressing the oscillations. 

Mathematically, weak damping means that all eigenvalues 
of A in (4) are still in the left half plane, whereas the principal 
eigenvalue which has the minimum absolute value has been 
very close to the imaginary axis. According to normal form 
analysis [25], the dynamics of system (4) can be decomposed 
into dynamics of individual modes, and we focus on the 
dynamics of the least stable mode described as below: 

i = (-7 -I- iu}o)z + (5) 

where 7 > 0 and 7 is close to zero. 

Let z = X + iy, ^ = ^x + i^y, and represent (5) in Cartesian 
coordinate, we have: 



The processes x{t), y{t) and z{t) are Ornstein-Uhlenbeck 
process which is stationary, Gaussian and Markovian. Lig la 
shows a typical sample path of x in system (6), where 7 = 
0.02, ujo = O.Stt, cr = 0.01. 

B. Limit Cycle 

Supercritical Hopf bifurcation, which leads to the emer¬ 
gence of a stable limit cycle, is another mechanism of power 
system oscillation. Unlike the weakly damped oscillation, the 
occurrence of Hopf bifurcation indicates that the equilibrium 
point of system has already lost its stability. The oscillation 
due to Hopf bifurcation can be regarded as an early warning of 
voltage collapse, since Hopf bifurcation usually occurs before 
saddle node bifurcation which results in the final voltage 
collapse [5]. Countermeasures like power reschedule, load tap 
changer blocking, and other emergency control may be needed 
to stop the voltage degradation. 

Limit cycles are inherently nonlinear phenomena that can 
be described only via nonlinear equations. Near the Hopf 
bifurcation point, the normal form of the stochastic system 
can be generally represented as: 

z = {y+ iu!h)z - \z\^z + (7) 

Here both z = x + iy coincides with the amplitude of the most 
unstable mode in the leading order and £, = + i^y is the 

same noise term as in linear case. It is possible to represent 
(7) in Cartesian coordinates: 



Lig. lb shows a typical sample path of x in system (8), where 
7 = 0.01, ujh = O.Stt, and a = 0.01. 
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C. Forced Oscillation 

Another type of power system oscillation is raised by 
forced oscillation from cyclic loads, control loops in power 
plants, diesel or hydro generators, wind turbines, etc. [13]- 
[15]. This kind of oscillation is not a result of the general 
dynamics of power system, instead, it is caused by an external 
forcing with a distinct oscillatory behavior. Particularly, forced 
oscillation may lead to voltage flickering when the frequency 
is around lOHz where the human eye is most sensitive; forced 
oscillation near the natural modes may result in resonance, and 
small disturbance is then amplified and expanded rapidly in 
the whole power system. The most effective countermeasure 
against forced oscillation is to locate and remove the external 
oscillation source. 

Mathematically, the system undergoing forced oscillation 
can be described as: 

i = (-7 -I- iujQ)z + (9) 

where wq is the natural frequency and fl is the forced fre¬ 
quency. 

Let z = x + iy, ^ and represent (9) in Cartesian 

coordinate, we have: 

f ^ \ ^ \ f ^ \ { Fcos{nt) +a^x '\ 

\y J -1 J \ y J V ^sin(m)-fcr^ J 

( 10 ) 

Fig. Ic shows a typical sample path of x in system (9) with 
7 = 1, Wo = 0.27r, F = 0.1, ft = O.Stt, and a = 0.01. 



timers) 



time(s) 


(a) 


(b) 



(c) 


Fig. 1: (a). A sample path of a; in a weakly damped system ( 6 ) 
with given parameters; (b). A sample path of x in system ( 8 ) 
with given parameters, which has a stable limit cycle; (c). A 
sample path of x in system ( 10 ) with given parameters, which 
has a forced oscillation. 


From Fig. la-lc, it is clearly seen that power system 
oscillations with different mechanisms all look similar to each 
other from time series data, which makes it hard to diagnose 
the exact cause and adopt the right control actions. In the 


next section, a novel method will be presented to identify 
the mechanisms for sustained oscillations by exploring the 
statistical signatures of time series data. 


III. STATISTICAL SIGNATURES OF DIFFERENT 
OSCILLATION MECHANISMS 


In this section, several statistical characteristics will be 
briefly introduced, and then elaborated analytically for differ¬ 
ent oscillation mechanisms. It will be shown that time series 
data of oscillation in different mechanisms exhibit different 
statistical signatures, based on which a diagnosis method is 
developed. 

Kurtosis is a descriptor of the shape of a probability 
distribution, and can be regarded as a measure of deviation 
from Gaussian distribution. Mathematically, kurtosis is defined 
as the standardized fourth moment around the mean of a 
distribution [26] [27]: 


E(a: — /i)^ /i 4 

(E(a: -/r)2)2 Var^[a:] 


( 11 ) 


where E is the expectation operator, y is the mean, ^14 is 
the fourth moment. The Gaussian distribution has a kurtosis 
of 3, so K = K — 3, termed as excess kurtosis, is often 
used so that the reference Gaussian distribution has a kurtosis 
of zero. In this paper, we use K in the subsequent analysis. 
For symmetric unimodal distributions, positive excess kurtosis 
indicates heavy tails and peakedness relative to the normal 
distribution, while negative excess kurtosis indicates light tails 
and flatness. In addition, the kurtosis is not affected by the 
variance since it is scaled with respect to the variance. 

The power spectral density (PSD) describes how the 
strength of a signal is distributed in the frequency domain. 
By Wiener-Khinchin theorem: 


SxH = 


p+oo 


Rx{r)t 


^ dr 


( 12 ) 


where Rx{t) = F(x{t)x{t + r)) is the autocorrelation func¬ 
tion. For example, periodic signals give peaks at a fundamental 
frequency and its harmonics, while white noise has a flat PSD. 

It will be shown that the kurtosis together with the PSD 
function can be used to distinguish the different mechanisms 
of oscillations. Firstly, the kurtosis can be used to distinguish 
the weakly damped oscillation from the limit cycle and the 
forced oscillation. Then, the PSD can be used to differentiate 
the limit cycle and the forced oscillation. Detailed analysis is 
to be presented for each scenario. 


A. Weakly Damped Oscillation 
System ( 6 ) can be represented as: 

q = + (13) 

where q = [cc, is an vector Ornstein-Uhlenbeck process, 

the PSD function can be calculated [28]: 

Sq{uj) = -;^{M + iw)~^ N N'^ {M'^ — iw)~^ 


( 14 ) 
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Specifically, the PSD of x is: 

a 2 ( 72+^2 


S'x(w) = 


+ UJ^) 


27r[(72 + Wq — oj'^Y + 47 ^ 0 ;^] 


(15) 



Fig. 2: Power spectral density of the sample path in Fig. la. 

Fig. 2 shows the PSD of the sample path presented in Fig. 
la. The power spectrum has a peak around /o = 0.15Hz with 
a wide bandwidth under the condition that 7 Ri 0. It indicates 
that if the system is weakly damped, the random excitation by 
the noise may lead to an oscillation near the natural frequency 
of the system. 

Regarding the kurtosis, the process x{t) is Gaussian since 
all linear functionals of a Gaussian process are also Gaussian. 
Hence, the kurtosis of x{t) should be approximately zero. 
Indeed, for the example shown in Fig. la, the kurtosis of the 
sample path is 0.24. 


B. Limit Cycle 

Consider system (8), near the Hopf bifurcating point, the 
amplitude of the limit cycle grows with y/j, and the angular 
frequency is approximately w/j. The solution can be approxi¬ 
mated as; 


{ x{t) 

V yW 


+ pi*)) cos{(j){t)) \ 


(16) 


where p = —2yp + ap^p is an Ornstein-Uhlenbeck process in¬ 
dependent of (j), = ojht + cTtp^tp is a Brownian motion with 

deterministic drift [29]. The PSD of x{t) can be represented 
as [30]: 

2 

5'x(w) = yF{al,u}h,u}) + ^F{al + 4y,ujh,u}) (17) 
where F is the PSD of cos{(j){t)): 


F{cr^,ujh,uj) 


al{u;^+ul+al/A) 

-Ujl-al/AY + 


4 2 


(18) 


The phase noise leads to a shift of the spectral peak, and even 
for very stable limit cycles (7 ^ tr^), the power spectrum 
has a non-vanishing bandwidth, leading to the temporal deco¬ 
herence which is known as “jitter” in electronic signal theory 
[31] [32]. Fig. 3 shows the PSD of the sample path in Fig. 
lb, which has a peak around 0.15Hz with a wide bandwidth. 

Even though the weakly damped oscillation and the limit 
cycle have similar power spectrums, they have different prop¬ 
erties with respect to the kurtosis. The kurtosis of x{t) = 



Fig. 3: Power spectral density of the sample path in Fig. lb. 


{■Xj + p{t)) cos{(j>{t)) is given below, refer to Appendix A for 
detailed derivation. 


Kurt[a;(f)] = 


3[(7-Var[p(f)])^-2(Var^[p(t)])] 
2(7 -I- Var[p(f )])2 


(19) 


From (19), we have that if Var[p(t)]) - k > {1 + V2), the 
kurtosis is negative; as k grows, i.e. the amplitude of the 
limit cycle is more and more larger than the noise intensity 
Var[p(t)], the kurtosis gets more and more closer to — which 
makes x{t) more distinguishable from Gaussian distribution, 
and thus from the weakly damped oscillation. For the example 
shown in Fig. lb, the kurtosis is -1.05, which is much farther 
away from zero compared with the weakly damped oscillation. 


C. Forced Oscillation 

We next consider system (10) whose solution can be repre¬ 
sented as: 

q = qo + qx (20) 


where go = [a^Ojt/o]^ satisfies 


q'o = Mqo + 


Fcos(Of) \ 
F sin(Hf) J 


( 21 ) 


and qi = [xi,yiY' is a vector Ornstein-Uhlenbeck process: 


q'l — Mqx + (22) 

We consider the real part x{t), the stationary solution (i.e., 
t —>■ -boo) of X can be represented as: x{t) = pcos(f2f) -|- 
xi{t), with the PSD function; 

Sxixi) = —p^[(5(a;o — H)-b i5(a;o + H)] (23) 

^ 2 ( 72 +^ 2 +^ 2 ) 

27r(72 -b Wq — -b 

Fig. 4 presents the PSD of the sample path shown in Fig. 
Ic. It demonstrates that in a well-damped system, the PSD 
function has a thin spike at the frequency ^ corresponding 
to the S function in (23). 

Compare Fig. 4 with Fig. 3 and Fig. 2, the forced oscillation 
is characterized by a thin spike in the power spectrum, which 
is distinct from the other oscillation mechanisms. In practice, 
the width of the spike will be determined by the decoherence 
time of the external forcing, so it will be narrow for truly 
periodic load variations. 
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Fig. 4: Power spectral density of the sample path in Fig. Ic. 


Furthermore, the kurtosis of x{t) = pcos(flf) + xi{t) can 
be represented as below, whose detailed derivation is given in 
Appendix A. 


Kurt [a: (f)] 


E[(x(f)-p)^] 
Var^ [x(t)] 


2 {l-\- o Var[aa(t)] ^2 
(24) 


(24) implies that the kurtosis of x{t) is negative, which 
means that the distribution of x{t) has light tails and 
flatness compared with Gaussian distribution. In addition, 
Kurt[a;(f)] > —and the absolute value of Kurt[a;(t)] 
depends on if the noise intensity of xi{t) is much 

less than the amplitude of the oscillation p, the kurtosis will 
be far away from zero and close to — which makes the 
distribution of x{t) distinguishable from Gaussian distribution. 
Indeed, for the example shown in Fig. Ic, the kurtosis is -1.44. 

Monte Carlo simulations are also done to demonstrate the 
characteristic of kurtosis for each mechanism. 500s simula¬ 
tion has been run for 100 times for each mechanism, and 
histograms of kurtosis value are shown in Fig. 5. The widely 
spread of kurtosis value makes different oscillation mecha¬ 
nisms readily distinguishable. Additionally, 90% confidence 
levels have been calculated. 90% confidence interval of the 
kurtosis for weakly damped oscillation is [-0.68, 0.56], 90% 
confidence interval of the kurtosis for limit cycle is [-1.29, 
-0.88], and that of the kurtosis for forced oscillation is [- 
1.45 -1.43]. The histograms and confidence intervals have 
demonstrated the pronounced distinguishability of kurtosis. 

From the above analysis, we see that different oscillation 
mechanisms exhibit different properties regarding the kurtosis 
and the PSD. Specifically, the kurtosis can be used to distin¬ 
guish the weakly damped oscillation, which has approximately 
zero kurtosis, from the forced oscillation and the limit cycle, 
which have negative kurtosis. Furthermore, the PSD can be 
used to distinguish the limit cycle which has a wide bandwidth, 
and the forced oscillation which has a thin spike. We hence 
propose an algorithm to diagnose the sustained oscillation 
causes as shown in Fig. 6. The threshold e has to be decided 
by considering noise intensity, confidence interval, system 
features and so forth in practical implementation; the detection 
of thin spike also depends on the tradeoff between noise 
intensity and forcing strength. Further studies are needed to 
improve the decision procedures in this algorithm. Besides, 
further investigations about control design are needed, while 
some preliminary results are to be presented in the following 
case studies, which is to locate the oscillation source or 



Fig. 5; Histograms of kurtosis value for all mechanisms. 



Fig. 6: Flowchart of the algorithm to diagnose sustained 
oscillations. 


problematic components via the kurtosis. 

IV. Case Studies 

Fig. 7 shows three different oscillation scenarios in 14-bus 
systems, from which the exact mechanism of each scenario can 
be hardly identified. In this section, we apply the proposed 
diagnostic method to identify the oscillation mechanism for 
each scenario using the simulated PMU data. Note that all the 
time series data is obtained from time-domain simulation of 
three different 14-bus systems which are modified from the 
IEEE 14-bus benchmark system. Stochastic loads are added 
to the systems. Particularly, all load fluctuations follow the 
Ornstein-Uhlenbeck process, with a = 0.01 in (3). Besides, 
we use the exponential recovery load which has been widely 
used in studying voltage stability [4] [33] and can be modelled 
as follows; 

P = kPoi^r (25) 
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Q = kQoi^f (26) 

>^0 

where fc is a dimensionless demand variable, Vq is the ref¬ 
erence voltage, and a and f3 are active and reactive power 
exponents depending on the type of load [33] [34]. All 
simulations were performed on PSAT-2.1.8 [34]. Parameter 
values of the base case are given in Appendix B, and all test 
systems are posted online; https://github.com/xiaozhew/Test- 
Systems-. 





500 1000 

time(s) 



(b) Voltage magnitude at Bus 2 


500 1000 

time(s) 


(c) Voltage magnitude at Bus 5 

Fig. 7: One specific voltage magnitude in three cases. 


A. Case I 

From the time series data shown in Fig. 7a, we see that the 
system is slowly evolving between 300s and 700s probably 
due to changing parameters. Given that the new steady state 
(between 700s to 1500s) suffers from oscillation, we intend to 
figure out the exact mechanism of the oscillation. 

Following the algorithm shown in Fig. 6, we firstly estimate 
the kurtosis of the time series data in Fig. 7a. The moving 
kurtosis of the time series data with a window size of 50s is 
presented in Fig. 8a. The peak in the figure corresponds to 
the transient dynamics due to changing parameters, while the 
kurtosis before and after the variation doesn’t change much. 
In fact, the kurtosis of the new steady state (between 700s 
to 1500s) is 0.006, which indicates that the random process 
is still approximately Gaussian and the nonlinearity doesn’t 
outstand. Therefore, we conclude that the sustained oscillation 
is due to weak damping. It implies that as the system evolves, 
the principal eigenvalue is getting closer to the imaginary axis, 
or in other words, the systems is getting closer to its stability 
boundary. 

We also examine the impact of measurement noise on the 
proposed oscillation diagnostic method. Measurement noise 
may wash out useful information and affect the observed 
statistics, and thus affect the detection and decision. To ad¬ 
dress this issue, the band-pass filter was applied to filter out 


Fig. 8: (a). The moving kurtosis of the voltage magnitude 
shown in Fig. 7a without measurement noise.; (b). The moving 
kurtosis of the voltage magnitude under the influence of 
measurement noise. 


measurement noise in [23]. It has been shown that the standard 
deviation (STD) of the measurement noise after Altering can 
be reduced to the order of 10“^ if the STD of the original 
measurement noise is 10“^. Following [23], we apply a white 
Gaussian noise with STD = 10“^ as the measurement noise 
to the simulated PMU data, and obtain the corresponding 
kurtosis shown in Fig. 8b. The estimated kurtosis of the new 
steady state (between 700s to 1500s) is -0.12. Even though the 
measurement noise deteriorates the statistic, it is still obvious 
from Fig. 8b that the kurtosis is close to zero, and thus system 
has weakly damped oscillation. 

The actual situation is that the variation of system voltages is 
a result of increasing loads. Starting from 300s, the exponential 
recovery loads at Bus 5 and Bus 12 increase gradually, and by 
380s, both loads grow by 8% and stop increasing afterwards. 
The increasing loads require more power support from the 
generators, making the field currents of generators closer to 
their limits. As the power required by the loads gets closer to 
the limit of the generators and the transmission network, the 
system is pushed closer to its voltage stability boundary, which 
leads to the weak damping of the complete power system. In 
order to maintain the stable operation of the system, damping 
of the system needs to be increased possibly by tuning the 
parameters of the excitation system. 

From this example, we see that the proposed oscillation 
diagnostic method is able to accurately distinguish the weakly 
damped oscillation from the other mechanisms, and provide 
important guidance to further design the right control actions. 

B. Case II 

In the second case, the time series data in Fig. 7b performs 
similarly as the previous case. Given that the new steady state 
exhibits oscillations after a transient period between 300s to 
700s, we want to diagnose the exact oscillation mechanism. 

By the proposed method, we firstly estimate the kurtosis 
to discern the weakly damped oscillation from the other 
mechanisms. Fig. 9a shows the moving kurtosis of the data in 
Fig. 7b, from which we can observe a significant change of the 
kurtosis before and after the variation. In fact, the kurtosis of 
the new steady-state time series data between 700s and 1500s 
is -0.58, which implies that the stochastic process has already 
deviated from Gaussian process. Hence, the weakly damped 
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oscillation has been eliminated from all the mechanisms in this 
scenario. Next, we estimate the PSD of the steady-state time 
series data between 700s and 1500s to differentiate the forced 
oscillation and the limit cycle. As presented in Fig. 9b, there 
is a clear peak around 0.12Hz in the power spectrum with 
a wide bandwidth. Hence, we conclude that the system has 
passed the Hopf Bifurcation point. The original equilibrium 
point has lost it stability and a stable limit cycle has emerged, 
which may be an early warning sign of voltage collapse. 



(a) (b) 


Fig. 9: (a). The moving kurtosis of the voltage magnitude 
shown in Fig. 7b; (b). The power spectral density of the steady- 
state time series data between 700s and 1500s in Fig. 7b. 

Similarly, we also incorporates measurement noise in this 
case. We apply a white Gaussian noise with STD = 10“^ 
as the measurement noise to the simulated PMU data, and 
obtain the corresponding moving kurtosis and PSD as shown 
in Fig. lOa-lOb. It is observed that both the kurtosis and 
the PSD perform similarly as the case without measurement 
noise. Particularly, the moving kurtosis of the new steady- 
state time series data between 700s and 1500s is -0.52, and 
the power spectrum has a clear peak around 0.12Hz with a 
wide bandwidth. Regarding the SNR, the amplitude of the 
limit cycle in the subspace of IV 2 I is 5 x 10“^ p.u., and 
variance of measurement noise is 10“®, hence the SNR is 
approximately —9.03 dB. In practice, both load fluctuations 
and electromechanical excitations act as an effective noise, 
thus making the SNR even lower. These effects are considered 
separately from the measurement noise in our study. 



Fig. 10: (a). The moving kurtosis of the voltage magnitude 
with measurement noise; (b). The power spectral density of 
the steady-state time series data between 700s and 1500s with 
measurement noise. 

In fact, the oscillation in this case is a result of continuously 
increasing the loads after 380s in Case I. The exponential 


recovery loads at Bus 5 and Bus 12 start to increase at 300s, 
grow by 12% at 420s, and stop increasing afterwards. The 
increasing loads result in the oscillation of over excitation 
limiter (OXL) for the generator at Bus 2, which further excites 
the fast variables of AVRs. The effect of OXL to voltage 
stability has been discussed in extensive literatures [4] [33] 
[35]- [38]. In this case, the competing effect between load 
dynamics and OXLs finally leads to the voltage instability. 
The power system passes the Hopf bifurcation point around 
400s when a stable limit cycle is born. 

To show the characteristics of kurtosis before and after Hopf 
bifurcation in this power system, Monte Carlo simulations 
have been done for both Case I and Case II. The histograms 
in Fig. 11 are plotted according to 55 360s samples of IV 2 I 
for both cases. The 90% confidence interval of kurtosis for 
Case I (i.e., weakly damped oscillation) is [-0.44, 0.10], and 
the 90% confidence interval of kurtosis for Case II (i.e., 
limit cycle) is [-0.82, -0.46]. The confidence intervals do not 
overlap with each other, which makes the decision procedure 
very easy. These results further validate the ability of kurtosis 
to distinguish between weakly damped oscillation and limit 
cycle. 



Fig. 11: Histograms of kurtosis value for all mechanisms. 

This example shows that the proposed oscillation diagnosis 
method is able to identify the limit cycle from all the mecha¬ 
nisms even when the SNR is low, and thus provides significant 
guidance to adopt the right control strategy in time. In this 
example, load increasing needs to be stopped to avoid further 
voltage degradation or even voltage collapse. 

We further observed that the kurtosis of real-time data 
may also help locate the oscillation source or the problematic 
components, because the absolute value of the kurtosis indi¬ 
cates the amplitude of oscillation under fixed noise intensity 
as shown in (19). By comparing the kurtosis of the voltage 
magnitudes at different buses which can be achieved from 
PMU data, the problematic components like the most stressed 
generator may be identified. In this example, the kurtosis of 
each voltage magnitude is shown in Table I. For illustration 
purpose, we do not consider measurement noise here. We see 










that the kurtosis of the voltage magnitude at Bus 2 has the 
maximum absolute value which indicates that the amplitude 
of voltage oscillation is larger at Bus 2 compared with other 
buses. The oscillation with larger amplitude implies more 
severe problem and potential oscillation source, which is true 
in this case. The activation of OXL of the generator at Bus 2 
results in the excitation of other fast variables of AVRs, and 
finally leads to the voltage oscillation issue. Additionally, as 
the power support for the load at Bus 5 is mainly from the 
generator at Bus 2, the load at Bus 5 needs to stop increasing 
in order to avoid further voltage degradation; parameters of 
the excitation system of the generator at Bus 2 also need to 
be tuned to suppress the oscillation. 


Bus 1 

-0.34 

Bus 8 

-0.22 

Bus 2 

- 0.54 

Bus 9 

-0.17 

Bus 3 

-0.05 

Bus 10 

-0.18 

Bus 4 

-0.23 

Bus 11 

-0.20 

Bus 5 

-0.22 

Bus 12 

-0.23 

Bus 6 

-0.19 

Bus 13 

-0.21 

Bus 7 

-0.23 

Bus 14 

-0.13 


TABLE I: The kurtosis of voltage magnitudes at different 
buses for the new steady state (between 700s and 1500s) in 
Case II. 


C. Case III 

In the last case, it is again hard to distinguish the oscillation 
mechanism directly from the time series data shown in Fig. 7c. 
Likewise, following the proposed method, we firstly estimate 
the kurtosis to discern the weakly damped oscillation from the 
forced oscillation and the limit cycle. The moving kurtosis of 
the time series data with a window size of 50s is presented 
in Fig. 12a, from which we can observe a significant change 
of the kurtosis before and after the variation. Particularly, the 
kurtosis of the new steady-state time series data (between 700s 
and 1500s) is -0.74, and hence the weakly damped oscillation 
has been ruled out from all the mechanisms. We next estimate 
the power spectral density of the steady-state time series data 
between 700s and 1500s in Fig. 7c. As presented in Fig. 12b, 
there is a thin spike around 0.15Hz in the power spectrum. 
Hence, we conclude that the system is experiencing forced 
oscillation at a frequency around 0.15Hz. 

We also include measurement noise in simulations. A white 
Gaussian noise with STD = 10“^ is added to the simulated 
PMU data, and the corresponding moving kurtosis and PSD 
are shown in Fig. 13a-13b. It is found that both the kurtosis and 
the PSD perform similarly as the case without measurement 
noise. Particularly, the moving kurtosis of the new steady-state 
time series data between 700s and 1500s is -0.56, and the 
power spectrum has a clear peak around 0.15Hz with a thin 
spike. Besides, the amplitude of the forced oscillation in the 
subspace of | V 51 is 2 x 10 “^ p.u., and variance of measurement 
noise is 10“®, hence the SNR is approximately —16.99 dB, 
which is very low and will be even lower if the influence of 
load fluctuations is considered. 



(a) (b) 


Fig. 12; (a). The moving kurtosis of the voltage magnitude 
shown in Fig. 7c; (b). The power spectral density of the steady- 
state time series data between 700s and 1500s in Fig. 7c. 



(a) (b) 


Fig. 13: (a). The moving kurtosis of the voltage magnitude 
with measurement noise; (b). The power spectral density of 
the steady-state time series data between 700s and 1500s with 
measurement noise. 


The actual situation is as follows, the exponential recovery 
load at Bus 5 increases by 5% between 300s and 350s, then 
one cyclic load joins at Bus 5 with a forced frequency 0.15Hz. 
Therefore, the increasing fluctuation of the bus voltage shown 
in Fig. 7c is due to forced oscillation. 

This example shows that the proposed diagnosis method 
can successfully identify the forced oscillation among all 
the mechanisms at a low SNR. As mentioned before. The 
mechanism diagnosis for the oscillation with small amplitude 
is necessary to implement mitigation control in a timely 
manner. 


Bus 1 

-0.33 

Bus 8 

-0.64 

Bus 2 

-0.48 

Bus 9 

-0.12 

Bus 3 

-0.29 

Bus 10 

-0.14 

Bus 4 

-0.60 

Bus 11 

-0.22 

Bus 5 

- 0.74 

Bus 12 

-0.30 

Bus 6 

-0.40 

Bus 13 

-0.27 

Bus 7 

-0.60 

Bus 14 

-0.12 


TABLE II: The kurtosis of voltage magnitudes at different 
buses for the new steady state (between 700s and 1500s) in 
Case III. 

Similarly, the kurtosis may help locate the oscillation 
source. The absolute value of the kurtosis indicates the ampli¬ 
tude of oscillation as shown in (24). By comparing the kurtosis 
of the voltage magnitudes at different buses, the location or 
region of the forced oscillation may be identified. In this exam- 
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pie, the kurtosis of each voltage magnitude is shown in Table 
II. For illustration purpose, we do not consider measurement 
noise here. We see that the kurtosis of the voltage magnitude 
at Bus 5 has the maximum absolute value, which indicates that 
the amplitude of voltage oscillation is larger at Bus 5 compared 
with other buses. It turns out that the cyclic load indeed is at 
Bus 5 with a forced frequency at 0.15Hz. Removing of the 
cyclic load at Bus 5 after the locating will then effectively 
eliminate the sustained oscillation. Further studies may be 
needed to design a more comprehensive method for locating 
the oscillation sources. 

V. VALIDATION ON DIFFERENT CASES 

In order to show that the proposed oscillation diagnostic 
method is robust, more simulation results are to be presented 
in this section. We change some parameters of the base case, 
and let the resulting system still experience Hopf bifurcation 
as the loads increase. Likewise, we also add a cyclic load 
to the base case in order to simulate forced oscillation. In 
particular, we set the active power of the exponential recovery 
load at Bus 13 to be 0.605 p.u. instead of 0.55 p.u. which 
is the value of the base case for Case /-///. In other words, 
we make the active power of the exponential recovery load 
increase by 10%. 




(a) Voltage magnitude at Bus 2 (b) Voltage magnitude at Bus 2 

1 . 051 -^^- 1 



1.0351 -'-^' 

0 500 1000 1500 

time(s) 

(c) Voltage magnitude at Bus 5 

Fig. 14; Voltage magnitudes in different mechanisms. 

Starting from 300s, the exponential recovery loads at Bus 
5 and Bus 12 increase at twice the speed of Case /, and by 
335s, both loads grow by 7% and stop increasing afterwards. 
Because the load increasing stop before the Hopf bifurcation, 
the system has weakly damped oscillation. Particularly, Fig. 
14a shows the voltage magnitude of Bus 2 of this case. 
Furthermore, if the increasing of exponential recovery loads 
at Bus 5 and Bus 12 do not stop until 360s, both loads grow 
by 12%. By then, the system has passed the Hopf bifurcating 
point, and the voltage magnitude at Bus 2 is shown in Fig. 14b. 
Regarding the forced oscillation, the manner of load changing 


is exactly the same as that in Case III, and Fig. 14c shows the 
voltage magnitude at Bus 5. 

We then follow the proposed oscillation diagnostic method 
and see whether the kurtosis and the PSD exhibit the expected 
characteristics. The moving kurtosis for the case before Hopf 
bifurcation is shown in Fig. 15a, from which it is observed 
that the kurtosis is still close to zero after the transient period. 
Actually, the kurtosis of the new steady state (between 700s 
to 1500s) is -0.24. As for the case that has passed the Hopf 
bifurcating point. Fig. 15b shows its moving kurtosis, from 
which it is seen that the kurtosis of the new steady state has 
been away from zero. In fact, the kurtosis of new steady state 
is -0.86. In the case of forced oscillation, the moving kurtosis 
is presented in Fig. 15c. It is seen that the kurtosis is away 
from zero, which is actually -0.51 for the new steady state. 
From these results, it is seen that kurtosis is a robust statistic to 
discern weakly damped oscillation from the other mechanisms. 



Fig. 15: (a). The moving kurtosis of the voltage magnitude 
shown in Fig. 14a; (b). The moving kurtosis of the voltage 
magnitude in Fig. 14b; (c). The moving kurtosis of the voltage 
magnitude shown in Fig. 14c; 


We further estimate the PSD for the cases of limit cycle and 
forced oscillation. The results are shown in Fig. 16a-16b. It is 
observed that the PSD of the limit cycle has a wide bandwidth, 
whereas the PSD of the forced oscillation exhibits a thin spike 
around 0.15HZ. Those results further validate that the PSD is 
able to distinguish the forced oscillation from the limit cycle. 

The extra simulation results given in this section are to 
demonstrate that the proposed oscillation diagnostic method 
is robust to the change of system parameters. 

VI. Conclusions and Perspectives 

This paper elaborates the mechanisms of power system 
oscillations in a unified mathematical framework, under which 
the statistical signatures of different oscillation mechanisms 
are investigated. Even though oscillations with different mech¬ 
anisms look alike in time series data, they exhibit distinct 
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(a) 


Fig. 16; (a). The power spectral 
of the time series data in Fig. 
density of the new steady state 
14c. 



(b) 


density of the new steady state 
14b. (b). The power spectral 
of the time series data in Fig. 


statistical signatures based on which an oscillation diagnosis 
method is developed. Particularly, the oscillation diagnosis 
method utilize the kurtosis to discern the weakly damped 
oscillation from the others, and then use the power spectral 
density to differentiate the limit cycle and the forced oscil¬ 
lation. Numerical example show that the proposed method 
can accurately identify the exact mechanism of sustained 
oscillation using PMU data, and also provide insights towards 
locating the oscillation sources. Though control actions are 
briefly discussed for each mechanism, further investigations 
are needed to better locate the problematic components and 
design the right control actions after identifying the exact 
oscillation mechanism. 


Appendix A 

Derivation of Kurtosis 

The kurtosis of x{t) = {-y/j+pit)) cos{(j){t)) can be derived 
as follows. Denote x(t) = y x(t)dt, then we have; 


E[(a;(f)'*] = sin'*(f2f) + sin^{i}t)xl{t) + xf(t) 


= gP +3p xf(t)+3(xf(t)) 
Therefore, the kurtosis of x(t) is; 

Kurt[a;(f)] = ~ - 3 = 

Var^[a;(f)] ^ (1 + 2^^)^ 


(33) 


(34) 


Appendix B 

PARAMETER VALUES OF THE BASE CASE 
The base case is modified based on the test system 
“d_014_dyn_mdr’ in PSAT-2.1.8 [34]. The line between Bus 
7-9 has been deleted. The active power of the load at Bus 2 
has been increased to 0.55 p.u. from 0.189 p.u., and the active 
power of the load at Bus 9 has been increase to 0.6 p.u. from 
0.413 p.u.. All loads are modelled as exponential recovery 
loads whose parameter values are given in Table III. Besides, 
there are two turbine governors for the generators at Bus 1 and 
2, with parameter values shown in Table IV. All generators are 
also controlled by OXLs whose parameter values are given in 
Table V. 


TABLE III; Exponential Recovery Load Parameter Values 


Parameter 

Value 

active power percentage kp 

ISo^ 

reactive power percentage kq 

loo^ 

active power time constant Tp 

lOs 

reactive power time constant Tq 

10 s 

static active power exponent as 

1 

dynamic active power exponent at 

1.5 for the load at Bus 4,5,9,11 

5 for the other loads 

static reactive power exponent Ps 

2 

dynamic reactive power exponent 0t 

2.5 for the load at Bus 4,5,9,11 
10 for the other loads 


TABLE IV; Turbine Governor Parameter Values 


E[a;(f)] = (y/j + p(t)) cos((f>(t)) = 0 (27) 

Var[a;(f)] = ((^7 -b p(f)) cos((/)(f)))2 = l(j + p2(t)) (28) 


E[(a;(f)'‘] = (72 -b 6jp‘^{t) +p‘^{t)) cos"‘((/)(f)) 

= + 67p2(b) + 3(p2(i))2) (29) 


Pai'ameter 

Value 

reference speed f 

1 p.u. 

droop R 

0.02 p.u. 

maximum turbine output 

1.2 p.u. 

minimum turbine output 

0.3 p.u. 

governor time constant Ts 

0 . 1 s 

servo time constant Tc 

0.45s 

transient gain time constant T 3 

Os 

power fraction time constant T 4 

12 s 

reheat time constant T 5 

50s 


The derivation of (29) uses the fact that xf (f) = 3(a;^(f))^ for 
Gaussian process and odd moments of Gaussian process are 
all zero. The kurtosis of x{t) is; 


Kurt [x{t)] 


3[h-pHt)r-2{p^m] 

2(7-bp2(i))2 


(30) 


Eollowing similar procedures for the limit cycle, the kurtosis 
of x{t) = pcos{ilt) + xiit) for the forced oscillation can be 
derived as follows. 


E[a:(<)] = pcos(nf)-b a;i(f) = 0 (31) 

_ n2 _ 

Var[a;(t)] = (pcos(r2f) + xi{t)y = — + xf(t) (32) 


TABLE V; Over Excitation Limiter Parameter Values 


Pai'ameter 

Value 

maximum held cuiTent 

for OXL 1-5 

5.1 p.u., 3.8 p.u., 2.5 p.u., 
3.2 p.u., 5.8 p.u. 

integrator time constant To 

30s for Generator 3 

10 s for the others 

maximum output signal 

100 p.u. 
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